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3D-Particle Tracking (3D-PTV) and Phase Sensitive Constant Temperature Ancmom- 
etry in pseudo-turbulence i.e., flow solely driven by rising bubbles — were performed to 
investigate bubble clustering and to obtain the mean bubble rise velocity, distributions 
of bubble velocities, and energy spectra at dilute gas concentrations (a ^ 2.2%). To 
characterize the clustering the pair correlation function G{r, 6) was calculated. The de- 
formablc bubbles with equivalent bubble diameter db = 4 — 5 mm were found to cluster 
within a radial distance of a few bubble radii with a preferred vertical orientation. This 
vertical alignment was present at both small and large scales. For small distances also 
some horizontal clustering was found. The large number of data-points and the non intru- 
siveness of PTV allowed to obtain well-converged Probability Density Functions (PDFs) 
of the bubble velocity. The PDFs had a non-Gaussian form for all velocity components 
and intermittency effects could be observed. The energy spectrum of the liquid velocity 
fluctuations decayed with a power law of —3.2, different from the ~ —5/3 found for 
homogeneous isotropic turbulence, but close to the prediction —3 by |Lance fc Bataille 



(19911 for pseudo-turbulence. 
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1. Introduction 



Bubbly pseudo-turbulence — i.e. a flow solely driven by rising bubbles — is relevant from 
an application point of view due to the omnipresence of bubble columns, e.g. in the 



chemical industry, in water treatment plants, and in the steel industry (Deckwer ( 1992 1) . 
A better understanding of the involved phenomena is necessary for scaling-up industrial 
devices and for optimization and performance prediction. This article wants to provide 
experimental data on pseudo-turbulence by means of novel experimental techniques. The 
main questions to be addressed are: 

(i) What is the preferential range and the orientation of bubble clustering in pseudo- 
turbulence? 

(ii) What is the mean bubble rise velocity and what kind of probability distribution 
does the bubble velocity have? 

(iii) What is the shape of the energy spectrum of pseudo-turbulence? 

1.1. Bubble clustering 

In dispersed flows, the hydrodynamic interaction between the two-phases and the particle 
inertia result in an inhomogeneous distribution of both particles and bubbles (see e.g. 
experimentally ( Ayyalasomayajula et a?.| [2006 Salazar et al. 2008; Saw et ar||2008| ) and 
numerically (Calzavarini et al. 2008a|& Bee et al. 2006a) and Toschi & Bodenschatz 



(20091 for a general recent review). Bubble clustering in pseudo-turbulence has been 



studied numerically (Smereka 1993 Sangani & Didwana 1993 Sugiyama et al. 2001 



Mazzitelli & Lohse 2009 ) and experimentally ( Cartellier & Riviere 2001 Risso & Ellingsen 



2002 Zenit et al. 2001 Roig & de Tournemine 2007[ ). The numerical work by Smereka 



(1993) and by Sangani & Didwana (19931 and theoretical work by van Wijngaarden 
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(19931, Kok (19931, and van Wijngaarden (20051 suggest that, when assuming potential 
flow, rising bubbles form horizontally aligned rafts. Three dimensional direct numerical 
simulations, which also solve the motion of the gas-liquid interface at the bubble's surface, 



have become available in the last few years. The work by Bunner & Tryggvason (2002a 



2003 1 suggests that deformability effect plays a crucial role for determining the orientation 



of the clustering. For spherical non-deformable bubbles these authors simulated up to 216 
bubbles with Reynolds numbers in the range of 12 — 30 and void fraction a up to 24%, 
and Weber number of about 1. For the deformable case, they simulated 27 bubbles with 
Reynolds number of 26, Weber number of 3.7, and a = 6%. The authors found that the 
orientation of bubble clusters is strongly influenced by the deformability of the bubbles: 
spherical pairs of bubbles have a high probability to align horizontally, forming rafts, 
whereas the non-spherical ones preferentially align in the vertical orientation. In a later 



investigation, where inertial effects were more pronounced, Esmaeeli & Tryggvason (2005 1 
studied both cases for bubble Reynolds number of order 100. In this case only a weak 
vertical cluster was observed in a swarm of 14 deformable bubbles. Their explanation 
for the weaker vertical clustering was that the wobbly bubble motion, enhanced by the 
larger Reynolds number, produces perturbations which do not allow the bubbles to align 
vertically and accelerate the break up in case some of them cluster. 

In spite of numerous experimental studies on pseudo-turbulence, bubble clustering 



has not yet been fully quantitatively analysed experimentally. Zenit et al. (20011 found 
a mild horizontal clustering using 2D imaging techniques for bubbles with particulate 



Reynolds number higher than 100. Cartellier & Riviere (2001 ) studied the microstructure 
of homogeneous bubbly flows for Reynolds number of order 10. They found a moderate 
horizontal accumulation using pair density measurements with two optical probes. Their 
results showed a higher probability of the pair density in the horizontal plane and a 
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reduced bubble density behind the wake of a test bubble. Risso & Ellingsen ( 2002 1 



performed experiments with a swarm of deformable bubbles ((4=2.5 mm), aspect ratio 
around 2, and Reynolds number of 800. They did not find clustering and suggested that 
in this low void fraction regime (a < 1.05%) there was a weak influence of hydrodynamic 
interaction between bubbles. 

1.2. Mean bubble rise velocity and statistics of bubble velocity 
The mean bubble rise velocity in bubbly flows is found to decrease with increasing bubble 



concentration whereas the normalized vertical fluctuation V Zirms /V z increases (Zenit 



et al. 2001 Bunner & Tryggvason 2003 Mercado et al. 2007). The interpretation is 



that when increasing the concentration the bubble-induced counterflow becomes more 
important and in addition the hydrodynamic interactions between bubbles become more 
frequent and hinder the upward movement of bubbles, provoking at the same time, an 
increment of the bubble velocity fluctuations. 

Next, Probability Density Functions (PDFs) of bubble velocities provide useful infor- 
mation for effective force correlations used in bubbly flow models in industry. PDFs in 



pseudo-turbulence have been obtained in the numerical studies of Bunner & Tryggvason 



(2002& 2003) and Esmaeeli & Tryggvason (20051. For non-deformable bubbles (Bunner 



& Tryggvason 2002&), the PDFs of velocity fluctuations have a Gaussian distribution. 



If deformability is considered ( Bunner & Tryggvason 2003 1 , the PDF of only one hori- 



zontal component of the velocity keeps a Gaussian shape while the non-Gaussianity in 
the PDFs was stronger at the lowest concentration a = 2%, recovering a Gaussian dis- 
tribution for a — 6%. However in that numerical work only a limited number of bubbles 
(Nf, = 27) could be considered and the different behavior in the two different horizon- 
tal directions reflect the lack of statistical convergence. Experimental PDFs of bubble 



velocities in pseudo-turbulence have been obtained by Zenit et al. (20011 and by Mer 
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cado et al. (20071. In these studies the bubble velocity was measured intrusively using 



an impedance probe technique. The amount of data-points used for the PDFs was not 
sufficient for good statistical convergence. 

1.3. Liquid energy spectrum 

In bubbly flows, one must distinguish between the energy input due to the bubbles and 
the energy input due to some external forcing — either of them can be predominant . [Lance| 



& Bataille ( 1991 ) suggested the ratio of the bubble-induced kinetic energy and the kinetic 



energy of the flow without bubbles as appropriate dimensionless number to characterize 



the type of the bubbly flow. Rensen et al. (2005) called this ratio the bubblance parameter 
b, defined as 



b = 



2 < 



(1.1) 



where a was already defined, U r is the rise bubble velocity in still water, and u' is 
the typical turbulent vertical fluctuation in the absence of bubbles. |Lance fc^ Bataille 



(1991) measured the liquid power spectrum in bubbly turbulence and observed a gradual 



change of the slope with increasing void fraction. At low concentrations the slope of 
the spectrum was close to the Kolmogorov's turbulent value of —5/3. By increasing a, 
the principal driving mechanism changed — the forcing now more and more originated 
from the bubbles and not from some external driving. In that regime the slope was 



close to —8/3. Having in mind equation (1.1 1 one may expect different scaling behavior, 



depending on the nature of the energy input that is dominant, namely b < 1 for dominant 
turbulent fluctuations or b > 1 for dominant bubble-induced fluctuations. Indeed from 
table 1 of Rensen et al. ( |2005[ ) one may get the conclusion that the slope is around —5/3 



for b < 1 and around — 8/3 for b > 1. Also Riboux et al. (2009) obtained a spectral lope 
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of about —3 in the wake of a swarm of rising bubbles in still water (6 = oo). Moreover, 



in numerical simulations Sugiyama et al. (2001 1 obtained the same spectral slope for the 



velocity fluctuations caused by up to 800 rising light particles, i.e. b = oo, with finite 



diameter (Re » 30). However, there are also counter-examples: e.g., Mudde et al. (19971, 



and Cui & Fan (20041 obtained around —5/3, though b = oo. Therefore, in this paper 
we want to re-examine the issue of the spectral slope in pseudo-turbulence (b = oo). 

1.4. Outline of paper 

The paper is structured as follows: in section [2] the experimental apparatus is explained 
and Particle Tracking Velocimetry technique, and Phase Sensitive Constant Temperature 
Anemometry are described. Section [3] is divided in three sub-sections: in the first part 
results on bubble clustering are shown, followed by the results on the mean bubble rise 
velocity and the bubble velocity distributions, and finally the power spectrum measure- 
ments are presented. Finally, a summary and an outlook on future work are presented in 
section |U 

2. Setup, tools, and methods 

2.1. Twente water channel 

The experimental apparatus consists of a vertical water tunnel with a 2 meter long 
measurement section with 0.45 x 0.45 m 2 cross section. A sketch depicting the setup 
is shown in figure [T] The measurement section has three glass walls for optical access 



and illumination (see |Rensen et al. (20051 for more details). The channel was filled with 
dcionized water until the top of the measurement section. The level of the liquid column 
was 3.8 m above the place where bubbles were injected. We used 3 capillary islands in 
the lowest part of the channel to generate bubbles. Each island contains 69 capillaries 
with an inner diameter of 500 /im. Different bubble concentrations were achieved varying 
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Figure 1: Experimental apparatus: a) Measurement section of length 2m, b) 4-camera 
PTV system, c) Capillary islands. 

the air flow through the capillary islands. We performed experiments with dilute bubbly 
flows with typical void fractions in the range 0.28% ^ a ^ 0.74% for PTV and in the 
range 0.20% ^ a ^ 2.2% for CTA. The void fraction a was determined using an U-tube 
manometer which measures the pressure difference between two points at different heights 



of the measurement section (see Rensen et al. 2005 ) . A mono-dispersed bubbly swarm 
with mean bubble diameter of 4— 5 mm was studied. Typical Reynolds numbers Re are of 
the order 10 3 , the Weber number We is in the range 2—3 (implying deformable bubbles) 
and Eotvos number Eo around 3 — 4. Table [I] defines these various dimensionless numbers 
and summarizes the experimental conditions. In figure [2] the mean bubble diameter d eq 
as a function of void fraction is shown. The values are within the range of 4-5 mm and 
show a slight increment with bubble concentration. 
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Figure 2: Equivalent bubble diameter d eq as a function of void fraction a. Standard de- 
viation as error bars. We measured the long di and short axis d s of 400 bubbles per 
concentration from 2D images. The equivalent diameter was obtained by assuming ellip- 
soidal bubbles with a volume equal to that of a spherical bubble, d eq = (dgdf) 1 ^ 3 



a(%) di(mm) V*(m/s) Re We Eo 
0.28 - 0.74 4- 5 0.16-0.22 1000 2-3 3-4 

Table 1: Typical experimental parameters: void fraction a, mean bubble diameter (db = 
\J dfd s ), di and d s are the long and short axis of 2D image of a given ellipsoidal bubble, 
mean rising velocity in still water V z , and Reynolds (Re = dbV g /uf), Weber (We = 
pdbV z /a), and Eotvos (Eo — pjgd\ja) numbers. Here Vf = 1 X 10~ 6 m /s is the 
kinematic viscosity of water, a — 0.072 N/m the surface tension at the air-water interface, 
and g the gravitational acceleration. 
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2.2. Particle Tracking Velocimetry 

In the last few years 3D-Particle Tracking Velocimetry (PTV) has become a powerful 
measurement technique in fluid mechanics. The rapid development of high-speed imaging 
has enabled a successful implementation of the technique in studies on turbulent motion 



Mordant et aZ.[2004 


Guala et aZ.|2005 


Bourgoin et al. 2006 


Berg et al. 



2006 Volk et al. 2008 Toschi & Bodenschatz 2009). The measured 3D spatial position 



of particles and time trajectories allow for a Lagrangian description which is the natural 
approach for transport mechanisms. 

Figure [l] also sketches the positions of the four high-speed cameras (Photron 1024-PCI) 
which were used to image the bubbly flow. The four cameras were viewing from one side 
of the water channel and were focused in its central region, at a height of 2.8 m above 
the capillary islands. Lenses with 50 mm focal length were attached to the cameras. We 
had a depth of field of 6 cm. The image sampling frequency was 1000 Hz using a camera 
resolution of 1024x 1024 pixel 2 . The cameras were triggered externally in order to achieve 
a fully synchronized acquisition. We used the PTV software developed at IfU-ETH for 
camera calibration and tracking algorithms. For a detailed description of this technique 



we refer to the work of Hoyer et al. (20051 and references therein. A 3D solid target was 
used for calibration. Bubbles were detected within a volume of 16x16x6 cm 3 with an 
accuracy of 400 fira. To illuminate the measuring volume homogeneous back-light and a 
diffusive plate were used in order to get the bubble's contour imaged as a dark shadow. 
The image sequence was binarized after subtracting a sequence-averaged background; 
then these images were used along with the PTV software to get the 3D positions of the 
bubbles. We acquired 6400 images per camera corresponding to 6.4 s of measurement 
(6.7 Gbyte image files). 

For higher bubble concentration, many bubbles are imaged as merged blobs and can not 
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Table 2: Number of detected bubbles and linked bubbles per image for all the concen- 
tration studied. For the highest void fractions there is a drop in the number of linked 
bubbles, since most bubbles are imaged as 2D merged blobs, which are not considered in 
the analysis. 



be identified as individual objects. These merged bubbles images are not considered for 
the analysis. Therefore the number Nb of clearly identified individual bubbles goes down 
with increasing void fraction. For the most dilute case (a = 0.28%) around w 190 
bubbles were detected in each image. This quantity dropped to nearly 100 for the most 
concentrated cases (a — 0.65% and a — 0.74% ). If a bubble is tracked in at least 3 
consecutive time-steps, we call it a linked bubble. Table [2] summarizes typical values of 
number of bubbles (Nb) and linked bubbles (iVjinfc) per image. 

2.3. Pair correlation function 

Particle clustering can be quantified using different mathematical tools like pair corre- 
lation functions (Bunner & Tryggvason 2002a), Lyapunov exponents ( Bee et aT||20066 1, 



Minkowski functionals (Calzavarini et al. 20086 1 , or PDFs of the distance of two con 
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Bubble i 



Figure 3: Definition of angle 9 between a pair of neighboring bubbles. 



sccutive bubbles in a time-series (Calzavarini et al. 2008a I. In this investigation the 



pair correlation function G(r, 8) is employed to understand how the bubbles are globally 
distributed. It is defined as follows: 



G(r, 6) 



V 



N b N b 



N b (N b - 1) 



(2.1) 



i=l j = 

where V is the size of the calibrated volume, Nb is the number of bubbles within that 
volume, Tij is the vector linking the centers of bubble i and bubble j, and r is a vector 
with magnitude r and orientation 9, defined as the angle between the vertical unit vector 
and the vector linking the centers of bubbles i and j , as shown in figure [3j From (2.1), the 
radial and angular pair probability functions can be derived. To obtain the radial pair 
probability distribution function G(r) one must integrate over spherical shells of radius 
r and width Ar, whereas for the angular pair probability distribution function G{9) an 
r-integration is performed. 
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2.4. Phase sensitive Constant Temperature Anemometry 

Hot-film measurements in bubbly flows impose considerable difficulty due to the fact that 
liquid and gas information is present in the signal. The challenge is to distinguish and 
classify the signal corresponding to each phase. The hot-film probe does not provide by 
itself means for a successful identification. Thus many parametric and non-parametric 
signal processing algorithms have been used to separate the information of both phases, 
i.e. thresholding ( Bruun|[l995 ) or pattern recognition methods ( Luther et q?.||2005 l. The 
output of these algorithms is a constructed indicator function which labels the liquid and 
gas parts of the signal. We follow an alternative way and measure the indicator function, 
therefore avoiding uncertainties when computing it. This can be achieved by combining 



optical fibers for phase detection to the hot-film sensor (see Cartellier k, Barrau 2001 
Julia et fl?.||2005) . 



The device, called Phase Sensitive CTA, was firstly developed by den Berg et al. (2009 1. 
In this technique, an optical fiber is attached close to the hot-film probe so that when 
a bubble impinges onto the sensor it also interacts with the optical fiber. Its working 
principle is based on the different index of refraction of gas and liquid. A light source is 
coupled into the fiber. A photodiode measures the intensity of the light that is reflected 
from the fiber tip via a fiber coupler. The incident light leaves the tip of the fiber when 
immersed in water, but it is reflected back into the fiber when exposed to air, implying 
a sudden rise in the optical signal. Thus the intensity of the reflected light indicates 
the presence of either gas or liquid at the fiber's tip. In this way, if both signals are 
acquired simultaneously, the bubble collisions can straightforwardly be detected without 
the need of any signal-processing method. For the construction of the probe we used 
a DANTEC cylindrical probe (55R11) and attached two optical fibers to it (Thorlabs 
NA=0.22 and 200 /im diameter core). We used two fibers in order to assure the detection 
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Figure 4: Phase Sensitive CTA probe with two attached optical fibers for bubble-probe 
detection. 

of all bubbles interacting with the probe. The fibers were glued and positioned at the 
side of the cylindrical hot-film, at a distance of about 1 mm from the hot-film. An 
illustration of the hot-film and fibers is shown in figure [4] The probe was put in the 
center of the measurement section positioning its supporting arm parallel to the vertical 
rising direction of the bubbles so that the axis of the optical probes are also aligned with 
the preferential direction of the flow, thus allowing for a better bubble-probe interaction 
and aiming at minimize the slow down of bubbles approaching the probe. |den Berg et al.\ 



(20091 measured flow velocity with and without the fiber being attached to the probe. 
They found that the presence of the fibers do not compromise the probe's bandwidth 
and that its influence on the power spectrum is negligible. 

With the signal of each optical fiber a discretized phase indicator function £ = 



can be constructed, whose definition follows Bruun (19951, namely 
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{1 liquid, 
(2.2) 
bubble. 

An unified phase indicator function is then obtained by multiplying the indicator 
functions of both fibers. A typical signal of Phase Sensitive CTA with two consecutive 
bubble-probe interactions is shown in figure |5a| There is one adjustable parameter to 
construct the indicator functions of the fibers. As explained above the phase discrimina- 
tion can be done by an intensity threshold of the optical fiber's signal. When the rising 
edge of the signal surpasses this threshold, the indicator function of the fiber must change 
from a value 1 to 0, as defined in equation |2. 2 1 If the optical fiber was positioned at the 
same place of the hot-film probe, the rising edge of its signal would occur precisely at the 
time when the bubble interacts with the hot-film. However, there is a separation of about 
1 mm between the fibers and the hot-film, so that the rising edge of the optical fibers' 
signal occurs actually with some delay. Therefore, to construct the phase indicator func- 
tion of the optical fibers, one must account for this delay and define the beginning of the 
interaction not where the signal surpasses the intensity threshold but some time before. 
The time used for this shifting was obtained considering the vertical separation between 
the optical fibers and the hot-film probe and mean bubble velocities: in our experiment 
a bubble travels 1 mm in about 5 to 7 ms. We shifted the beginning of the collision 8 
ms prior to the optical fibers' signal starts to rise from its base value. The shifting value 
was double checked by analyzing the histograms of the duration of bubble collisions. As 
it can be observed from figure |5b[ with this shift sometimes part of the CTA signal when 
the bubbles is approaching the probe (< 10 ms) was lost, but we noticed no effect on the 
spectrum when varying the shift duration, provided the bubble spikes were still removed. 

With the phase indicator function only pieces containing liquid information of the 
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Figure 5: Typical phase Sensitive CTA signal already transformed into velocity, (a) two 
consecutive bubble-probe interactions for a=0.2%, (b) zoom in of a bubble collision. 

time series are used for further analysis — i.e. the part of the signal between the two 
bubble interactions shown in figure |5a| For each part of the time-series containing liquid 
information the spectrum was calculated. Then all spectra were averaged and the power 
spectrum for each bubble concentration was obtained. In this way neither an interpolation 
nor auto-regressive models for gapped time series were necessary. In our experiments the 
gas volume fractions varied from 0.2 ^ a ^ 2.2%. The CTA was calibrated using a 
DANTEC LDA system (57N20 BSA). The calibration curve was obtained by fitting the 
data points to a 4th order polynomial. The total measuring time was 1 hour for each 
concentration and the sampling rate of the hot-film and optical fibers was 10 kHz. 
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3. Bubble clustering and distributions in pseudo-turbulence 

3.1. Radial pair correlation 

The radial pair correlation G(r) was obtained for all bubble concentrations studied. 
Figure [6] shows G(r) as a function of the normalized radius r* = r/a, where a is a mean 
bubble radius. As shown in figure [2] the mean equivalent bubble diameter is within the 
range 4-5 mm. Therefore we can normalize with one mean bubble radius and we picked 
a = 2 mm for all concentrations. We observe in figure [6] that the highest probability to 
find a pair of bubbles lies in the range of few bubble radii r* ss 4 for all concentrations. 
The probability G{r) of finding a pair of bubbles within this range increases slightly 
with increasing a. For values r* < 2 one would expect that G(r) = 0. However, in our 
experiments we found G(r) ^ for r* < 2, due to the fact that the bubbles are ellipsoidal 
and deform and wobble when rising. 

We now estimate the error bar in the correlation function G(V), originating from incom- 
plete bubble detection, as seen from table [2j With increasing a the fraction of detected 
bubbles decreases. For a = 0.28% we detect Nb w 200 so one would expect for a=0.74% 
to detect N b w (0.74/0.28) x 200 w 500 but we are detecting only 110, «20% of them. 
In order to investigate the reliability of the pair correlation function due to this loss of 
bubble detection, we studied a set of randomly distributed particles. We generated 500 
particles at 6000 times and calculated the radial and pair correlation functions, obtaining 
G(r) = 1 and G(9) = 1 as expected for randomly distributed particles. Then we kept 
only 20% of the particles and recalculated the correlation functions to see how much 
they deviate from 1. We found that the maximum deviation was less than 0.1 for both 
the radial and the angular correlation functions. This is much smaller than the structure 
revealed in the G(r) curve in figure [6j We therefore consider the clustering with the pre- 
ferred distance of r* = 4 as a robust feature of our data, in spite of incomplete bubble 
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Figure 6: Radial pair probability G(r) as a function of dimensionless radius r* . The 
different curves are for different bubble volume concentrations a (given in percent). 

detection. In figure [6] this error bar corresponding to a maximum error of G(r) at the 
most concentrated flow is also shown. 

3.2. Angular pair correlation 

The orientation of the bubble clustering was studied by means of the angular pair corre- 
lation G(8) using different radii for the spherical sector over which neighboring bubbles 
are counted. Figures [7a] [7b] and [7c] show the results for r*=40, 15, and 5, respectively. 
The plots were normalized so that the area under the curve is unity. For all radii and 
concentrations, pairs of bubbles cluster in the vertical direction, as one can see from the 
highest peaks at 8/ir = and 9 /it = 1. The value of 9/ir = means that the reference 



18 J. Martinez, D. Chehata, D.P.M. van Gils, C. Sun, and D. Lohse 

bubble (at which the spherical sector is centered) rises below the pairing one. For 9/tt = 1 

the reference bubble rises above the pairing bubble (see figure |3| . When decreasing the 
radius of the spherical sector, i.e. when probing the short range interactions between the 
bubbles, we observe that a peak of the angular probability near tt/2 starts to develop. 
The enhanced probability at this angle range is even more pronounced in figure [7c| where 
the peak of G(8) for horizontally aligned bubbles is just slightly lower than that for ver- 
tical clustering. It is worthwhile to point out that the vertical alignment of the bubbles is 
very robust and is present from very large to small scales, as the angular correlation for 
different spherical sectors is always higher at values 9/ir=0 and 1 than at value 0/ir=0.5. 
The maximum error bar for the angular correlation for the most concentrated flow at 
a = 0.74%, when 20% of particles are detected, was 0.1 as explained above and is also 
shown in figure [7j It is much smaller than the structure of the signal. 



For comparison, we consider again the work of Bunner & Tryggvason (20031, who 
found that pairs of bubbles have a higher probability to align vertically, though for a 



much higher concentration (a = 6%) than employed here. Bunner & Tryggvason (20031 
found that the vertical alignment was not as robust as in our case, since with increas- 
ing r* the angular correlation at and n became less dominant. Another significant 
difference between the findings of our experimental work and their simulations is that 
horizontal alignment was more pronounced with larger radii of the spherical sector and 
not when decreasing r* . Our experimental results clearly show the main drawback and 
today's limitation when solving the flow at the particle's interface: the simulations are 
still restricted to a small number of particles, which is not sufficient to reveal long range 
correlations. 
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Figure 7: Normalized angular pair probability G{6) as a function of angular position 9/ir 
for various bubble concentrations a (see insets) and three different bubble-pair distances: 
(a) r* = 40, (b) r* = 15, and (c) r* = 5. Maximum error bar for the angular correlation 
of 0.1 at a=0.74%. 
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3.3. Interpretation of the clustering 

What is the physical explanation for a preferred vertical alignment of pairs of bubbles 
in pseudo-turbulence? Through potential flow theory, the mutual attraction of rising 



bubbles can be predicted (Batchelor 19671, the application of potential theory to our 



experiments remains questionable (van Wijngaarden 19931, as we are in a statistically 
stationary situation where bubbles have already created vorticity. Our findings are con- 
sistent with the idea that deformability effects and the inversion of the lift force acting 



on the bubbles are closely related to the clustering. Mazzitelli et al. ( 2003 1 showed nu- 
merically that it was mainly the lift force acting on point-like bubbles that makes them 
drift to the downflow side of a vortex in the bubble wakc|f] Furthermore, when account- 
ing for surface phenomena, Ervin & Tryggvason ( 1997[ ) showed that the sign of the lift 
force inverts for the case of deformable bubbles in shear flow so that a trailing bubble is 
pulled into the wake of a heading bubble rather than expelled from it. In such a manner 
vertical rafts can be formed. Experimentally some evidence of the lift force inversion has 



been observed by Tomiyama et al. ( 2002 1 as lateral migration of bubbles under Poiseuille 
and Couette flow changed once the bubble size has become large enough. Numerical 
simulations of swarm of deformable bubbles without any flow predicted a vertical align- 



ment (Bunner & Tryggvason 2003). An alternative interpretation of the results, due to 
Shu Takagi (private communication (2009)) goes as follows: small, pointwise, spherical 
bubbles have a small wake, allowing for the application of potential flow. The bubbles 
then horizontally attract, leading to horizontal clustering. In contrast, large bubbles with 
their pronounced wake entrain neighboring bubbles in their wake due to the smaller pres- 
sure present in those flow regions, leading to vertical clustering. Further efforts are needed 
to identify and confirm the main mechanism i.e., cither lift or pressure reduction in the 



t See figure 2 in 



Mazzitelli et al. 



(20031 sketching the dynamics. 
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bubble wake — leading to a preferential vertical alignment, for example through experi- 



ments with small, spherical, non-deformable bubbles as achieved by Takagi et al. (20081 
through surfactants or with buoyant spherical particles. 

3.4. Average bubble rise velocity 
Bubble velocities were calculated by tracking the bubble positions which were linked in 
at least three consecutive images. The mean bubble rise velocity can thus be obtained 
as a function of the bubble concentration. Figure [8] shows the three components of the 
bubble velocity; the coordinate system corresponds to the one depicted in figure [l] The 
terminal rise velocity of a single bubble with d(, = 3.9 mm and with the same water- 
impurity condition is also shown in figure [8] it has a value of 0.26 m/s. A decrease in the 
mean bubble rise velocity with concentration is observed in our experiments within the 
experimental error showed in figure § The mean bubble rise velocity is 0.22 m/s for the 
most dilute case a = 0.28% and decreases until a value of 0.16 m/s for a = 0.51%. 
The interpretation of this finding is that in this parameter regime the velocity-reducing 



effect of the bubble-induced counterflow (see |Batchelor 1972 ) and the scattering effect 



overwhelm the velocity-enhancing blob-effectjf] ( |Brenner 1999[ ) , implying that a blob of 



rising bubbles rises faster than a single one. For values of a ^ 0.56% the mean values are 

again larger, around 0.18 m/s. This increment of the mean rise velocity could be a result 

of our experimental error. As mentioned in section |2.2[ the number of detected bubbles 

at higher concentrations decreases by a factor 3 as compared to the most dilute cases. 

To check whether this increment was coming from detection of blobs of bubbles rather 

than single ones, we did experiments with a single camera positioned perpendicularly 

to the flow. The 2D images were used to track bubbles manually making sure that 

the trajectories indeed corresponded to single bubbles. In figure [8] the mean bubble rise 
f Originally suggested for sedimenting particles 
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Figure 8: (a) Mean bubble velocity as a function of bubble concentration a. All three 
components are shown: o V x ; □ V y ; T mean bubble rise velocity V z ; A mean bubble 
rise velocity obtained by particle tracking from single camera 2D images; ♦ terminal rise 
velocity for a single bubble with d(, = 3.9 mm. The error bars were obtained by estimating 
the 95% confidence interval for the mean. 

velocity from the one-camera 2D analysis are plotted. A similar behavior is observed, first 
a decrease with concentration, followed by a slight increase for the most concentrated 
flows, confirming the 3D PTV results. 

3.5. Bubble velocity distributions 

Figure [9] shows the PDF of all velocity components at the most dilute concentration 
(a = 0.28%). The number of data-points used for calculating the PDFs was larger than 
9x 10 4 for the highest concentration a=0.7% and of order 3x 10 5 for the most dilute case. 
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Even for the most concentrated flow, the number of data-points was large enough to 

assure the statistically convergence of the PDFs. All PDFs show non-Gaussian features, 



as nicely revealed in the semi-log plot of the PDFs (figure 10 1. 

In order to quantify the non-Gaussianity of the PDFs, the flatnesses of the distributions 
were calculated. Their values are shown in the inset of figure [9] and are within the range 
6—13. The flatness of the vertical component has always the highest values in the whole 
range of void fraction. 



The effect of the concentration on the PDFs is also shown in figure 10 As in that 
figure, the PDFs are shown on a semi-logarithmic scale so that the deviation from the 
Gaussian-like shape become more visible, it is revealed that there is no substantial change 
in the shape of the distributions with increasing bubble concentration. 

We would like to compare the non-Gaussianity of the PDFs ( "intermittency" ) found 
here with a comparable system, namely Rayleigh-Benard convection as the analogy be- 



tween bouyancy driven bubbly flows and free thermal convection is interesting (see Cli- 



ment & Magnaudet 19991. In Rayleigh-Benard convection a fluid between two parallel 



plates is heated from below and cooled from above, see Ahlers et al. (20091 for a recent 
review. Prominent coherent structures in this system are thermal plumes, which arc fluid 
particles either hotter or colder than the background fluid. Due to the density difference 
with the background fluid, hotter plumes ascend and colder ones descend. The system is 
solely buoyancy driven. Particularly, for large Pr the plumes keep their integrity thanks 
to the small thermal diffusivity, so that the similarity with pseudo-turbulence is appeal- 

Does the statistics of the velocity fluctuations in Rayleigh-Benard share a similar 
f We stress however that of course there are differences between the two systems which have 



been pointed out by Climent & Magnaudet 



(1999) 
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Figure 9: Probability Density Functions for the three bubble velocity components for 
a = 0.28% normalized with the rms values of each component. V x : dotted line with 
squares, V y : dashed line with circles (both with V mean velocities as expected), V z : dotted- 
dashed line with triangles. The respective black solid lines (without symbols) show a 
Gaussian with the same mean and width as the measured distributions. The inset presents 
the flatness of the distribution as a function of the concentration a. The horizontal solid 
line in the inset represents the flatness for a normal distribution. Squares, triangles, and 
circles have the same meaning as in the main plot. 

behavior with that of bubbles in pseudo-turbulence? In Rayleigh-Bcnard convection, 
the PDF of the vertical velocity fluctuations of the background fluid — i.e. the central 



region of the cell — exhibits a Gaussian distribution (Daya & Ecke 2001 1 . Qiu & Tong 



On bubble clustering and energy spectra in pseudo-turbulence 



25 



Q 
D- 



10° 



10- 



10" : 



1 ci- 



io-' 



10" f 



-15 



-10 







I 

■ 0.28% 
X 0.35% 






▼ 0.41% 
O 0.51% - 
A 0.56% 
♦ 0.65% 
+ 0.74% 
gaussian - 


»*• / 

XT / 


\ <* 

*** 

ll O I 



-5 5 

v x' v x,rms 



10 



15 



Q 
D- 



10° 



10"' 



10" : 



10" : 



10- 



10" ; 



-15 



Q 
D. 



10° 



10- 



10" ; 



10" : 



10"' 



10" 1 



1 

(b) / 




1 

■ 0.28% 
X 0.35% 
T 0.41% 

O 0.51% - 
A 0.56% 
♦ 0.65% 
+ 0.74% 
Gaussian- 


^ / 

/ 

>« ♦ / 





-10 



10 



15 



VyVy.rms 



I 




i 


I 

■ 0.28% 


(c) 






X 0.35% 
▼ 0.41% 
0.51% - 






£ 4 


A 0.56% 
♦ 0.65% 

+ 0.74% 
gaussian- 




/ 








i i \ 





-15 



-10 



10 



15 



v z /v z , rms 



Figure 10: Same Probability Density Functions of the bubble velocity normalized with 
the rms values of each component as in figure [9] but now on a semi logarithmic scale to 
better reveal intcrmittency effects for various concentrations. As a reference, a Gaussian 
curve with the same mean and standard deviation as for the distribution with a = 0.28% 
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(20011, and Sun & Xia (20051 measured the vertical velocity fluctuations in the region 
where the plumes abound, they found that the PDF still follows a Gaussian function. 
Those measurements, carried out for Pr sa 4, indicate clearly that the PDF of the plume 
velocity fluctuations can be described by a Gaussian function, which differs significantly 
from the statistics of the bubble velocity in pseudo-turbulence. A possible reason of this 
difference could be that buoyancy in pseudo-turbulence is much stronger than that of 
the plumes in Rayleigh-Benard system. 

3.6. Energy spectra of pseudo-turbulence 

Figure [TT] shows the energy spectra for all gas fractions. As it can be seen, the slope of 
the energy spectrum hardly depends on the volume fraction — all curves show a slope of 
about —3.2. We stress that this scaling behavior is maintained for nearly 2 decades, much 



wider than it had been reported in prior observations ( Lance fc Bataille||1991 Sugiyama 
et q?.|200T Riboux et a?.|2009 1 of this steep slope of pseudo-turbulence spectra. As it was 
mentioned in section |2~4| the way the power spectrum was calculated in this investigation 
differs from previous ones in two aspects: First, the indicator function has been measured 
by means of the optical fibers and second, an energy spectrum has been calculated for 
all individual liquid segment, before the final spectrum is obtained through averaging. 

One wonders whether the duration of our interrupted time series is large enough to 
resolve the low frequency part of the spectrum: if the duration of these segments is too 
short, then indeed the low frequencies in the power spectrum can not be resolved. On the 
other hand, if the duration of the liquid segments is large enough, then all frequencies 
in the spectrum are well resolved. Figure [12] shows the distribution of the logarithm of 
the non-interrupted time series duration for three different concentrations. For a=2.2% 
(the most concentrated bubbly flow with more bubble-probe interaction, thus shorter 
liquid segments) around 90% of the segments used to construct the spectrum have a 
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Figure 11: Liquid energy spectra for different void fractions. • a— 0.2%, ▲ a=0A%, ■ 
a=0.8%, ♦ a=1.3%, T a— 1.7%, + a=2.2%>. The arrow shows the onset frequency of the 
scaling. 

duration larger than 0.05 s. Comparing with figure |11[ we can see that for frequencies 
higher than 20 Hz we have resolved the inertial range, where the slope is « —3.2. Thus 
the measurement of the spectra is consistent. 

Why does the slope differ from the Kolmogorov value —5/3? One might expect a dif- 
ferent scaling in pseudo-turbulence, as the velocity fluctuations are caused by the rising 
bubbles and not by the Kolmogorov- Richardson energy cascade, initiated by some large 
scale motion. The difference between these two scalings is not yet completely understood, 
but there are some hypothesis on its origin. One possible explanation for the different 



scaling in pseudo-turbulence was given by Lance & Bataille (19911. They argued that 



eddies from the bubble wake are immediately dissipated before decaying towards smaller 
eddies, which would lead to the —5/3 scaling. They derived a —3 scaling, balancing the 
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Figure 12: Distribution of the time duration of the parts containing liquid information in 
the CTA-signal used to calculate the spectrum for three different bubble concentrations. 
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The solid vertical line corresponds to the onset frequency of the scaling in figure 11 
arrow in that plot. 

spectral energy, and assuming that the characteristic time of spectral energy transfer is 
larger than those of dissipation and production. More evidence that wake phenomena are 



related to the —3 scaling has been given by Risso et at. (2008) and by Roig & de Tourne 



mine 



(20071. They showed that bubbles' wakes decay faster in pseudo-turbulence than 



in the standard turbulent case with the same energy and integral length scale. They 
also proposed a spatial and temporal decomposition of the fluctuations in order to gain 
more insight into the different mechanisms. In their experiments they used a fixed array 



of spheres. Very recently Riboux et al. (20091 measured the spatial energy spectrum in 



pseudo-turbulence by means of PIV. They measured the spectrum just immediately af- 
ter a bubble swarm has passed and obtained a —3 decay for wavelengths larger than the 
bubble diameter (2 mm < l c <7 mm). For wavelengths smaller than the bubble diameter 
they found that the spectrum recovered the —5/3 scaling. Their findings showed that 
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the scaling is independent on the bubble diameter and void fraction in their range of 

parameters investigated (0.2% ^ a ^ 12% and db = 1.6 — 2.5 mm), which we also find. 
Their conclusion is that the —3 scaling is only a result of the hydrodynamic interactions 
between the flow disturbances induced by the bubbles. Their arguments for this state- 
ment are that the scaling appears for wavelengths larger than the bubble diameter and 
that a different scaling was found for smaller wavelengths. It is worthwhile to emphasize 
that in our case, we measured within the swarm where production is still maintained and 
steady. The —3 scaling in our measurements is in the range of 8 mm down to hundreds of 
micrometer^f] thus for lengths not only larger than db but also for those up to one order 
of magnitude smaller than it. This supports that the dissipation of the bubble wake is 



involved and still a valid explanation for the —3 scaling as proposed by Lance fc Bataille 



(19911. In any case, the papers by Sugiyama et al. (20011; Roig & de Tournemine (20071; 



Risso et al. (20081; Riboux et al. (20091 and the present one show that the —3 scaling is 



typical for pseudo-turbulence. 

One further result supporting this idea is the one obtained by [Mazzitell i fc Lohse| 



( 2009 1 who performed numerical simulations of microbubbles in pseudo-turbulence mod- 
eling them as point particles. Their DNS treated up to 288000 bubbles, where near- 
field interactions were neglected, thus wake mechanisms can not be accounted for, and 
effective-force models were used for the drag and lift forces. They obtained a slope of the 
power spectrum close to —5/3 typical for the turbulent case. This gives evidence that 
the bubble's wake missing in the point particle approach — and its dissipation play a 

very important role for the —3 scaling of the energy spectrum in pseudo-turbulence, 
f estimated by considering the mean bubble rise velocity and the starting and ending fre- 



quencies of the scaling of the spectrum in figure |TT| 
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4. Summary 

We performed experiments on bubble clustering using 3D-PTV. This is the first time 
that the technique is used to investigate bubbly flows in pseudo-turbulence at very dilute 
regimes (a < 1%). Bubble positions were determined to study bubble clustering and 
alignment. For that purpose the pair correlation function G(r, 6) was calculated. As the 
radial correlation G(r) shows, pairs of bubbles cluster within few bubble radii 2.5 < 
r* < 4. Varying the bubble concentration does not have any effect on the clustering 
distance. The angular pair correlation G(8) shows that a robust vertical alignment is 
present at both small and large scales, as it is observed when varying the radius of the 
spherical sector (r*=40, 15, and 5). Decreasing the radius of the spherical sector shows 
that horizontal clustering also occurs, as the peak of the angular correlation around it/2 
starts to grow with decreasing values of r* . 

Probability density functions of bubble velocity show that all components of bub- 
ble velocity behave differently from Gaussian. The implementation of this non-intrusive 
imaging technique assures enough data points to obtain convergence in PDFs. The im- 
provement achieved in the number of data-points, compared with previous experimental 
investigations, is of order 10 2 . Furthermore, this allowed us to show the intermittent sig- 
nature that bubble distributions have for all components. The flatness values for these 
velocity distributions are in the range of 6—13. The distribution of the rise velocity 
showed the highest values of flatness, e.g. « 13 at a = 0.28%. The non-Gaussianity can 
be a result of the cluster formation mechanism, where the rise velocity of single bubbles 
are affected by the faster collective motion of clusters. However further investigations are 
needed to fully understand its origin. 

We have shown that the power spectrum in pseudo-turbulence (b = oo) decays expo- 



nentially with a slope near —3 which is consistent with the theoretical scaling that Lance 



On bubble clustering and energy spectra in pseudo-turbulence 



31 



k, Bataille ( 1991 1 derived and supporting the hypothesis that bubble wake mechanisms 



are closely related to it. We have shown that the implementation of Phase Sensitive CTA 
for studying bubbly flows is of great advantage, allowing for a direct recognition and 
discard of bubble-probe interactions. 

The next step of our research will be to investigate bubble clustering for b << 1, 
where turbulent effects become dominant. Another line of research is to analyze pseudo- 



turbulence with smaller bubbles (e.g. achieved by surfactants (Takagi et al. 2008)), to 
study the effect of the length scale of the bubble on the spectra. 
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